function [assigned, matched] = func_DA(type, pref_stu, pref_sch, capacities)

% DA   Implements Gale-Shapley Deferred Acceptance algorithm
%
% Inputs:
% - type: 
%   'stu': student-proposing DA (student optimal)
%   'sch': school-proposing DA (school optimal)   
% - pref_stu: = transforming all utilities to be positive; students' `ranking' of schools (J x I matrix sorted by school_id and student_id). Highest rank = most preferred / 0 if unranked.
% - pref_sch: school's ranking of students (J x I matrix sorted by school_id and student_id). Highest rank = most preferred / 0 if unranked.
% - capacities: school capacities 
%
% Output:
% - matched: JxI matrix of school_id x student_id matches
% - assigned: each student's assigned school_id (I x 1 sorted by student_id)

% alternative 1 faster if sample if small / alternative 2 faster for large sample (e.g. asymptotic)
alternative = 2;

% -------------------- %
% FOR TESTING FUNCTION %
% -------------------- %

% Set test = 1 to test function for student proposing DA / 2 for school proposing DA
test = 0;

%%% Change this
if test == 3
    type='stu';
    pref_stu = Unconstrained.Stu_rk_pref(:,:,1)';
    pref_sch = Unconstrained.Ranks(:,:,1);
    %ranking_sch = fliplr(Ranks2);
    capacities = Capacities;
end
if test == 1
    type='stu';
    pref_stu = pref_stu2;
    pref_sch = Ranks;
    ranking_sch = fliplr(Ranks2);
    capacities = Sch_capacities;
end
if test == 2
    type='stu';
    pref_stu = pref_stu2;
    pref_sch = Ranks;
    ranking_sch = fliplr(Ranks2);
    capacities = sch_capacities_as;
end
 
% ---------- %
% PARAMETERS %
% ---------- %

% # of students, # of schools
[nsch, nstu] = size(pref_stu);

% ----------------------- %
% SCHOOL OPTIMAL MATCHING %
% ----------------------- %

if strcmp(type,'sch') == 1
    
    % Logical array to select last admitted student in each school
    offers_selector =  zeros(nsch, nstu); 
    for ss=1:nsch
        offers_selector(ss,capacities(ss))=1;   
    end
    offers_selector = logical(offers_selector);
    
    stop = 0;
    while stop == 0; 
        % Sort students by priority in each school
        [priorities_stu, ~] = sort(pref_sch,2,'descend');
        priorities_stu = priorities_stu';
        % Priority of last admitted student
        prior_last_admitted_stu = priorities_stu(offers_selector');
        % Offers are made to all students with higher priority
        offers = (pref_sch>=repmat(prior_last_admitted_stu,[1 nstu]));
        % Add existing match to the matrix
        % Each student chooses the school with the highest utility out of offers
        matching.sch = (pref_stu.*offers == repmat(max(pref_stu.*offers),[nsch 1])).*offers;
        % Declined offers
        declined_offers = offers - matching.sch;
        % Stopping criterion: new declined offers
        stop = 1-any(declined_offers(:));
        % Update the preferences of schools that have been declined by students
        pref_sch(logical(declined_offers)) = 0;
    end

end

% ----------------------- %
% STUDENT OPTIMAL MATCHING %
% ----------------------- %

if strcmp(type,'stu') == 1
    
    % First: create JxI matrix where (j,i) is the id of the student that school i ranks in i-th position (in descending order). col1 = most preferred students
    [~, ranking_sch] = sort(pref_sch,2,'descend');    
    
    % Logical array to select last admitted student in each school
    offers_selector =  zeros(nsch, nstu); 
    for ss=1:nsch
        offers_selector(ss,capacities(ss))=1;   
    end
    offers_selector = logical(offers_selector);

    stop = 0;
    while stop == 0;        
        % School to which students propose: most preferred among non-rejected
        proposals = double(pref_stu ~= 0 & pref_stu == repmat(max(pref_stu),[nsch 1]));
        % Schools evaluate proposals
        if alternative == 1
            % Rank of all proposing students (col no = stud id)
            stu_ranked = pref_sch.*proposals;
            [evaluate, ~] = sort(stu_ranked,2,'descend');
            evaluate = evaluate';
            % Priority of last admitted student
            cutoff = evaluate(offers_selector');
        elseif alternative == 2
            % Loop over schools
            for ss=1:nsch
                % Id of proposing students
                id_proposing = ranking_sch(ss,:).*proposals(ss,ranking_sch(ss,:));
                % remove non-proposing students
                id_proposing(id_proposing == 0) = [];
                % Priority of last admitted student
                % If remaining seats
                if size(id_proposing,2)<capacities(ss)
                   cutoff(ss,1) = 0;                
                % If excess demand (or exactly filled)
                elseif size(id_proposing,2)>=capacities(ss)
                    % iD of last admitted student
                    id_last_admitted = id_proposing(capacities(ss));
                    % Cutoff = priority of last admitted student
                    cutoff(ss,1) = pref_sch(ss,id_last_admitted);
                end
            end
        end        
        % Any student below cutoff is rejected
        rejected_stu = proposals.*(pref_sch<repmat(cutoff,[1 nstu]));
        %disp(mat2str(sum(sum(rejected))))
        matching.sch = proposals - rejected_stu;
        % Update the preferences of rejected students
        pref_stu(logical(rejected_stu)) = 0;
        stop = 1-any(rejected_stu(:));
    end
end

% ---------------- %
% MATCHING OUTCOME %
% ---------------- %

% matched: JxI matrix of school x student matches
matched = matching.sch;

% assigned: 1xI vector of studen'ts assigned school (0 if unassigned)
assigned = (max(matching.sch.*repmat((1:nsch)',[1 nstu]),[],1))';
             